Anti-ferromagnetic ordering in arrays of superconducting 7r-rings 
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We report experiments in which one dimensional (ID) and two dimensional (2D) arrays of 
YBa2Cu307-i-Nb 7r-rings are cooled through the superconducting transition temperature of the 
Nb in various magnetic fields. These 7r-rings have degenerate ground states with either clockwise 
or counter-clockwise spontaneous circulating supercurrents. The final flux state of each ring in the 
arrays was determined using scanning SQUID microscopy. In the ID arrays, fabricated as a single 
junction with facets alternating between alignment parallel to a [100] axis of the YBCO and rotated 
90° to that axis, half-fluxon Josephson vortices order strongly into an arrangement with alternating 
signs of their magnetic flux. We demonstrate that this ordering is driven by phase coupling and 
model the cooling process with a numerical solution of the Sine-Gordon equation. The 2D ring 
arrays couple to each other through the magnetic flux generated by the spontaneous supercurrents. 
Using 7r-rings for the 2D flux coupling experiments eliminates one source of disorder seen in similar 
experiments using conventional superconducting rings, since 7r-rings have doubly degenerate ground 
states in the absence of an applied field. Although anti-ferromagnetic ordering occurs, with larger 
negative bond orders than previously reported for arrays of conventional rings, long-range order is 
never observed, even in geometries without geometric frustration. This may be due to dynamical 
effects. Monte-Carlo simulations of the 2D array cooling process are presented and compared with 
experiment. 



C 

O 
o 



> 

CO 
CN 
CO 

o 
in 
o 

-I— > 



i 

C 

o 
o 



X 



I. INTRODUCTION 



In superconducting rings, the requirement of a single- 
valued, macroscopic quantum mechanical wave-function, 
combined with the intimate relation between the quan- 
tum mechanical phase and the vector potential, results 
in flux quantization. 1 * 2 - Under the appropriate conditions, 
two of the flux quantized states can become degener- 
ate, characterized by time-reversed, persistent supercur- 
rents circulating in a clockwise and counter-clockwise di- 
rection. These macroscopic circulating currents can be 
thought of as an analog to an electronic spin. In pio- 
neering work, Davidovic et ali^ showed that arrays of 
superconducting rings can be used as models for spin 
systems, since they interact anti-ferromagnetically upon 
cooling. However, the Davidovic arrays never showed 
long-range anti-ferromagnetic ordering. They speculated 
that one reason for this lack of ordering is that, for there 
to be two degenerate states, such rings must be cooled in 
a field equivalent to a half-integer multiple of the super- 
conducting flux quantum $o = h/2e per ring. Therefore 
variations in the lithographically patterned areas of these 
rings result in different fluxes through different rings in 
the same field, causing disorder. Superconducting 7r-rings 
have an intrinsic phase shift of 7r in the absence of an 
externally applied field or supercurrent^*§*£ Such rings 
have a doubly degenerate ground state in zero applied 
flux, should not have this source of disorder, and may 
therefore be a more ideal model for a spin system. 

However, use of 7r-rings for a model spin system re- 
quires a large number of rings. The first supercon- 
ducting 7r-ringSf2i2ii2iii which depended on the momen- 
tum dependence of the pairing wavefunction in the 



high-T c cuprate perovskite superconductors, were made 
using technologies which would be difficult to extend 
to many devices. Recently a technology that allows 
for photolithographic fabrication of 7r-shift devices and 
arrays of great complexity, using YBa2Cu307_5-Au- 
Nb (YBCO-Au-Nb) ramp-edge tunneling junctions has 
been demonstrated. 12,13 « Moreover, it has recently been 
demonstrated that 7r-rings can also be fabricated using 
Josephson junctions with ferromagnetic layers in the tun- 
nel barriersii^i^i^ In this paper we report on experi- 
ments in which the YBCO ramp-edge technology was 
used to fabricate large arrays of 7r-rings. A first report 
on work with similar arrays appeared in Ref . Il7l The ar- 
rays were cooled in various magnetic fields, and the final 
"spin" states of the arrays were determined with scan- 
ning SQUID microscopy. Long range anti-ferromagnetic 
(AFM) ordering was observed in the ID arrays. Although 
stronger anti-ferromagnetic correlations were observed in 
our 2D 7r-ring arrays than were reported previously for 
the conventional (0-ring) arrays, in neither case did AFM 
ordering extend beyond a few lattice sites. Although us- 
ing 7r-rings eliminated one source of disorder, inhomo- 
geneous flux biasing due to the applied fields required 
for 0-rings, there are other sources of disorder, including 
inhomogeneities in the ring critical currents and criti- 
cal temperatures. We will discuss the influences of these 
sources on our results, as well as dynamic effects, using 
Monte-Carlo modelling of the cooling process. 



II. EXPERIMENTAL METHODS 

Various samples consisting of one- and two- 
dimensional 7r-ring arrays have been realized using 
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FIG. 1: Schematics of the 3 types of facetted junctions stud- 
ied in this work, (a) A continuous, single facetted junction, 
(b) Junction with the individual half-fluxons electrically iso- 
lated from one another, (c) Double row of continuous facetted 
junctions. 



ramp-type YBa2Cu307_«5 - Au - Nb Josephson contacts. 
The fabrication of ramp-type YBa2Cu307_5 - Au - 
Nb Jose phson has been described previously in detail 
in Ref.'s Il2ll3l In short, the samples were prepared 
by first epitaxially growing a bilayer of [001]-oriented 
YBa 2 Cu307_(5 and SrTi03 by pulsed-laser deposition on 
[001]-oriented SrTiC>3 single crystal substrates. For the 
ID array samples a 150 nm YBa2Cu307_,5 and a 100 nm 
SrTi03 film were used, while for the 2D array samples 
a 340 nm YBa2Cu307_,5 and a 67 nm SrTi03 film were 
used. In these bilayers the basic layout of the structures, 
which will be described below in more detail, is defined 
by photolithography and Ar ion milling. This process 
results in interfaces with a slope of 15 — 35° with respect 
to the substrate, which provides access to the afe-planes 
of the YBa 2 Cu307_5 and allows the exploitation of 
d-wave phase effects. Special care is taken to align all 
interfaces accurately along one of the YBa2Cu307_,5 
[100] axes. After etching the ramp and cleaning of the 
sample, a 7 nm YBa2Cu307_5 interlayer is deposited, 
the function and properties of which are described in 
Ref. [13, followed by the in-situ pulsed-laser deposition 
of a Au barrier-layer of 12 nm for the ID array samples 
and 20 nm for the 2D array samples. A 160 nm Nb 
counter electrode is then formed by dc sputtering and 
structured by lift-off. Each chip contained several 
reference-junctions. At T = 4.2 K, these showed a 
typical critical current per micrometer junction- width 
I c /w ~ 0.10 mA//im for the ID array samples and 
I c /w w 0.18 mA//im for the 2D array samples. From 
these, a value for the Josephson penetration depth 
Xj ps 1/im (T = 4.2 K) for all samples is deduced, which 
is the characteristic length scale over which Josephson 
vortices (fluxons) extend. 

The sample magnetic fields were imaged with a high 
resolution scanning SQUID microscope i^i^&SL The 
SQUID microscope images shown here were made at a 
temperature of T < 5K, with the sample cooled and 
imaged in the same fields. The size of the pickup loop 
used will be indicated for each image. The samples were 



FIG. 2: Schematic diagram of two types of 7r-rings used for 
the 2D arrays reported in this paper. Type (a) was used in 
the square arrays, and type (b) was used in the triangular, 
honeycomb, and kagome arrays. 
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FIG. 3: Scanning electron microscopy images of 4 arrays of n- 
rings, with 2.7^im junctions and 11.5^im ring to ring spacings. 

warmed through the superconducting critical tempera- 
ture of the Nb and cooled at controlled rates either using 
a non-inductive heater, or by passing warm He gas past 
the sample. 

A first configuration for which the generation and 
coupling of half-integer flux quanta was investigated is 
the zigzag array^ several instances of which are shown 
schematically in Figure ^ I n these structures, the d- 
wave order parameter of the YBa2Cu307_a induces a 
difference of tt in the Josephson phase-shift <fi across the 
YBa2Cu307_5-Au-Nb barrier for neighboring facets. For 
facet lengths a in the wide limit, i.e., a >> Xj, the lowest- 
energy ground state of the system is expected to be char- 
acterized by a spontaneous generation of a half-integer 
flux-quantum at each corner. This half-fluxon provides 
a further 7r-phase change between neighboring facets, ei- 
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ther adding or subtracting to the d-wave induced 7r-phase 
shift, depending on the half-flux quantum polarity. In 
both cases this leads to a lowering of the Josephson cou- 
pling energy across the barrier, as this energy is propor- 
tional to (1 - coscfi ). We studied three types of facetted 
junctions. The first (Fig. QJi) was an isolated, continu- 
ous junction with many adjacent facets. The second (Fig. 
njj) was lithographically patterned to electrically isolate 
each half-fluxon from its neighbor. The final type (Fig. 
|TJ0 had two continuous facetted junctions close together, 
but electrically isolated from one another, to test for field 
coupling between ID arrays of half-fluxons. 

The 2-dimensional 7r-ring arrays were made up of indi- 
vidual rings patterned as indicated schematically in Fig- 
ure|3 There were two types of rings, square (Fig. Et) and 
hexagonal (Fig. Eb)- The rings were patterned into ar- 
rays with 4 different geometries, as shown in the scanning 
electron microscopy images of Figure |3| In all cases the 
nearest- neighbor distances in these arrays was 11.5/xm 
center to center. The details of the ring geometries, crit- 
ical currents, self-inductances, and mutual-inductances 
are given in Table 1. 



III. RESULTS 

Figure 01 shows representative scanning SQUID mi- 
croscope images of 6 zigzag — tt facet junctions. All 
of these images were of samples on the same substrate, 
and imaged in the same cooldown in nominal zero field ( 
< 0.5/iT) at T = 4.2K. In all three types of facetted junc- 
tions, half-fluxon Josephson vortices were spontaneously 
generated at the points where the facets met, as expected. 
When the connected ID arrays (Fig. [TJi) were cooled in 
zero field through the niobium superconducting transi- 



TABLE I: Details of 2D lattice samples; critical current den- 
sity is J c « 5 x 10 8 A/m 2 





Square 


Triangular 


Honeycomb 


Kagome 


height YBCO (nm) 


340 


340 


340 


340 


height STO (nm) 


67 


67 


67 


67 


height Nb (nm) 


160 


160 


160 


160 


width JJ 1 {lira) 


2.75 


2.70 


2.70 


2.70 


width J J 2 (/im) 1 


2.75 


3.75 


3.75 


3.75 


hole (/im) 


2.50 


~3.10 


~3.10 


~3.10 


self-inductance (pH) 2 


3.93 


~3.90 


~3.90 


~3.90 


self-inductance (pH) 3 


4.57 








nearest neighbor 
distance (/im) 


11.5 


11.5 


11.5 


11.5 


nearest neighbor 
mutual (pH) 3 


0.025 









1 designed value 

2 estimate with standard formulas 

3 estimate using FastHenry, Ay_sco = 160nm and Ajvf,=40nm 
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FIG. 4: Scanning SQUID microscope images of zigzag YBCO- 
Au-Nb — 7r junctions, cooled in nominally zero field, and im- 
aged with a 4pm diameter pickup loop, (a) Continuous junc- 
tion with 40^im facet lengths, (b) Electrically disconnected 
junction with 40^m between facet corners, (c) Electrically 
disconnected junction with 20/im between facet corners. Two 
parallel facetted — n junctions with 40/xm (d), 20/xm (e), 
and 10/xm (f) between facet corners. The apparent curvature 
of the junctions in these images is an artifact of the scanning 
mechanism. 



tion temperature, the signs of the half-fluxons strongly 
tended towards perfect anti-ferromagnetic ordering, in 
which the persistent supercurrents alternated between 
clock- wise and counter clock- wise flow (Fig. 0^,d,e,f). 
There are two possible mechanisms for this ordering. In 
the first, ordering is driven by the minimization of the to- 
tal junction Josephson coupling energy during the cool- 
ing process. We refer to this as "phase" coupling, because 
the Josephson currents and energies are determined by 
the difference in superconducting phase across the junc- 
tion (j). In the second, the anti-ferromagnetic ordering 
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FIG. 5: Scanning SQUID microscope images of a facetted 
YBCO-Nb — 7T junction with 10 facets each 40/im long, 
cooled in fields of OnT (a), 32nT (b), 74nT (c), and HOnT 
(d), and imaged at 4.8K with an Sfim square pickup loop. 



is driven by a minimization of the total magnetic field 
energy in the junction and its environment. We refer to 
this as "field" coupling. In order to determine the rel- 
ative strengths of these mechanisms, facetted junctions 
were fabricated with each half-fluxon electrically isolated 
from its neighbor (Fig. H 3 ). This should eliminate the 
phase coupling mechanism. Indeed, when this is done the 
anti-ferromagnetic ordering in these ID arrays is much 
weaker. Two examples are shown in Fig. 0Jd,c. In Fig. 
0Jd the half-fluxons all have the same orientation. We be- 
lieve that this is because they align with a small residual 
field. Fig. shows a more random alignment of the 
isolated half-fluxons. 

In order to test for magnetic field coupling between the 
ID half-fluxon chains, double facetted junctions (Fig. [TJ;) 
were also fabricated and tested. It appeared that also in 
this case the phase coupling was stronger than the field 
coupling: Fig. EJl shows a section of a 40/im facet dou- 
ble junction that shows in-phase alignment between the 
two anti-ferromagnetically ordered ID chains. Because 
this arrangement places the positive half-fluxons in the 
lower left chain closest to the negative half-fluxons in the 
upper right chain, this is the lowest energy arrangement. 
However, in sections of the ID chains which show defects, 
as in the center of the upper right chains in Fig. and 
Fig. 0f, the interchain alignment goes from in-phase to 
out-of-phase when the interchain ordering has a defect, 




FIG. 6: SQUID microscopy images of 4 electrically discon- 
nected arrays of 7r-rings, in the geometries illustrated in Fig. 
El with values for the ring parameters given in Table I. These 
images were taken at 4.2K with a 4fim diameter pickup loop 
after cooling in nominally zero field at 1-lOmK/sec. The 
spin up fractions x+ and full scale variations in the scanned 
SQUID sensor flux (A$ s ) were £+=0.56, A$ s = 0.12$ for 
the square lattice; £+=0.51, A$ s = 0.11 ( I > o for the honey- 
comb lattice; £+=0.50, A<3> 3 = 0.10<3?o for the triangular lat- 
tice; and £+=0.54, A$ s = 0.12<I>o for the kagome lattice re- 
spectively. 



but one chain does not develop a second defect to align 
the interchain spins. Therefore it appears that the en- 
ergy cost to create a defect is larger than the energy gain 
from making neighboring chains in-phase. 

When the facetted junctions are cooled in an exter- 
nally applied magnetic field, one spin direction becomes 
energetically favored over the other, but there is a com- 
petition between this energy and the anti-ferromagnetic 
coupling energy during the cooling process. An example 
is shown in Figure [SJ More detailed results and mod- 
elling of this facetted junction will be described in the 
next section. 

In the electrically disconnected 2D lattices of Fig. [3] 
the square and honeycomb arrays are geometrically un- 
frustrated, as their magnetic moments can be arranged 
so that all nearest neighbors have opposite spins. In 
contrast, the triangle and kagome lattices are geomet- 
rically frustrated, since it is impossible for all of the 
rings to have all nearest neighbors anti-ferromagnetically 
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FIG. 7: SQUID microscopy images of a honeycomb array in 
the geometry illustrated in Fig. |3J cooled in nominally zero 
field through the Nb superconducting transition temperature 
at different cooling rates. Each panel is labelled with the 
cooling rate and experimentally determined bond order a. 
The white circles superimposed on the images label 6-ring 
loops in the honeycomb lattice in which the rings are perfectly 
anti-ferromagnetically ordered. 



100/im 




FIG. 8: Difference image obtained from subtracting the image 
of Fig. |7Ji from Fig. |7|d, to determine which rings flipped sign 
after successive cooldowns. 



aligned. Fig. shows examples of scanning SQUID mi- 
croscope images of the arrays of Fig. [3] after cooling in 
nominally zero field. Although regions of perfect anti- 
ferromagnetic ordering are seen in the unfrustrated ar- 
rays (Fig. E|i,b), anti-ferromagnetic ordering beyond a 
few lattice distances was never observed. Nevertheless, 
antiferromagnetic correlations were seen in all the 2D tt- 
ring arrays. A measure of the short range antiferromag- 
netic correlations is the bond order 



XAF 

2x+X- 



(1) 



where xaf is the fraction of the nearest neighbor pairs 
with opposite supercurrent circulation, and x + (x_) is the 
fraction of rings which have up (down) moments. Perfect 
antiferromagnetic correlation would correspond to a = - 
1. The minimum possible bond order at zero applied field 
for the frustrated triangular and kagome lattices is a = 
-1/3^ The images in Fig. Elare labelled with values for 
the bond orders. 

It is to be expected that the anti-ferromagnetic order- 
ing of the 2D arrays should improve if they are cooled 
more slowly through the Nb superconducting transition. 
Figure shows SQUID microscope images of the same 
region of the honeycomb 2D lattice of Fig. |3 after cool- 
ing at various rates. The individual panels are labelled 
with the cooling rates and final state bond orders. The 
anti-ferromagnetic ordering increases with slower cooling 
rates. One question that can be asked is: Do particular 
regions of the 2D array order more strongly than oth- 
ers? The white circles in Fig. outline the 6-member 
rings in the honeycomb arrays in which all neighbors are 
anti-ferromagnetically aligned. This provides a conve- 
nient way of visualizing regions of local order. It appears 
that there are no correlations between the positions of 
the ordered 6-member rings from cooldown to cooldown, 
and we conclude that the ordered regions are randomly 
distributed in space. 

In the 0-ring experiments of Davidovic et ali&i re- 
peated cooling resulted in a particular ring often being in 
the same final state (spin-up or spin-down). This is pre- 
sumably the result of the rings having slightly different 
effective areas, and therefore cooling in slightly different 
effective fields, since 0-rings must be cooled in finite fields 
for there to be degenerate states. The individual rings in 
our 2D 7r-ring arrays appear to have random final states. 
Figure |H1 shows a difference image between Fig. [7Jd and 
Fig. UJ:, which were taken after successive cooldowns. 
In this image, the rings which did not switch sign from 
one cooldown to the next are not visible, the rings which 
switched from down to up appear black, while the rings 
which switched from up to down appear white. Of the 
roughly 663 rings in this field of view, 152 switched from 
down to up, and 162 switched from up to down. This is 
consistent with random switching during cooldown. 
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FIG. 9: Modelling of the cooldown of the facetted junction of 
Fig. |S] The junction has 10 faces, each 40 microns long, with 
alternating 0- and 7r-intrinsic phase shifts. 
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FIG. 10: Results from cooldown of the facetted junction of 
Fig. |^| in several externally applied magnetic fields. Half- 
fluxons with magnetic fields in one direction are represented 
by open symbols, those with the opposite sense with closed 
symbols. The left two panels are experimental results from 
two successive runs; the right panel represents modelling as 
described in the text. 



IV. MODELLING 

For the current purposes we treat the facetted ramp 
edge junction as a linear junction with alternating regions 
of 0- and ir- intrinsic phase shifts 0(x) extending in the re- 
direction, with the junction normal in the z direction, and 
the junction width w in the y direction small compared 
with the Josephson penetration depth Xj = ^/h/2e^odjc, 
where d is the spacing between the superconducting faces 
making up the junction, and j c is the Josephson critical 
current per area of the junction. The quantum mechani- 
cal phase drop cf>(x) across the junction is the solution of 
the Sine-Gordon equation: 



^ 1 • (Af \ 



9(x)). 



(2) 



Analytical solutions to this equation are availabler 2 * 2 ^ 2 ^ 
but here we solve Eq. [3numericallyi 25 i 26 i 27 i 28 i 29 Defining 
a dimensionless coupling parameter 



A 2 



(Ax) 2 



(3) 



the differential equation Eq. [21 turns into a difference 
equation on a grid of size Ax: 



1 



o 



(4) 



Taking the net current across the junction equal to zero, 
using a reduced externally applied magnetic field h e = 



2ed\jH e /fi, and using the relation between the gradient 
of the phase and the field H in the junction: 



(5) 



we find boundary conditions that are described as differ- 
ence equations, where rij is the total number of junctions: 



(6) 



These coupled difference equations are solved using a re- 
laxation method to find the solution 4>{x), iterating to 
convergence. As the junction cools through the super- 
conducting transition temperature T c , the supercurrent 
density j c increases from zero, so that Xj decreases from 
infinity. To model the cooling process, we first solve Eq. 
2] for A j much larger than the facet length Lf, then de- 
crease Xj by a small amount, take the previous solution 
as the starting point for the next solution, and repeat un- 
til A j << Lf. An example is shown in Fig. |5| Here the 
externally applied flux was set at $ e = H e dL = $o/2, the 
total junction length _L=400/im, Ax=l/im, with 10 facets 
each of length L/=40/im, an initial Aj = 100/im, which 
was decreased to 4/zm in 24 equal steps. The field thread- 
ing the junction is given by H(x) = ($>o/2nd)d4>/dx. The 
final solution strongly favors anti-ferromagnetic arrange- 
ment of the half-fluxon "spins" . 

We can understand why the — 7r facetted junctions 
cool into perfect anti-ferromagnetic order, while the iso- 
lated linear arrays of 7r-rings do not, by considering the 
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FIG. 11: Calculated energies for the honeycomb ring arrays 
of Fig. |3 as a function of the reduced temperature. Eb is 
the barrier to thermally activated flipping of the sign of the 
spontaneous magnetization, E c is the strength of the spin-spin 
coupling energy, and ksT is the thermal energy. 



energetics of the cooling process. Written as a difference 
equation, the free energy for a particular state of the 
facetted junction becomes: 



F v = 



hj c wAx 
2e 



71 j —1 

£ 

i 



- COs((j) n + 0„) + - (4>n+l ~ 4>n) 2 

(7) 

Numerical solution of the Sine-Gordon equation as de- 
scribed above for the facetted junction of Fig. [5] shows 
that the free energy/facet, when the junction is in per- 
fect anti-ferromagnetic ordering, is ~ -5x 10 5 K/Aj(^m). 
The energy cost to form a defect, by flipping one spin, 
is - 7.2xl0 5 e- 45/A/( ' im) if. This implies that when the 
free energy /facet is comparable to k^T, the energy cost 
to form a defect is «7.2xl0 5 K : it is energetically favor- 
able to form perfect anti-ferromagnetic ordering in linear 
— 7r facetted junctions cooled in zero field. 

Figure ITU1 compares the results from repeated cooling 
of the 10-facet zigzag junction under various magnetic 
fields with modelling using the numerical solution of the 
sine-Gordon equation outlined above. Note that there are 
some disorder effects in the cooling process, as evidenced 
by the slight differences between the two experimental 
runs. The qualitative features of the data are reproduced 
by the modelling, although the experimental results are 
not as symmetric with respect to inversion in position, or 
with respect to field reversal, as predicted. One possible 
source of the observed asymmetry might be field gradi- 
ents. However, putting a linear field gradient into the 
model did not improve the fit with experiment. Another 
source of asymmetry might come from the asymmetry of 
the junction and lead geometry. 

The cooling process for isolated 7r-rings is different 



than for electrically connected 7r-rings. For simplicity 
we model our isolated 7r-rings considering a single induc- 
tance L and Josephson junction critical current I c . The 
details of the formulas will be different for a two-junction 
ring, 30 but we do not expect the physics to be qualita- 



tively different, 
be written as 



The free energy U of such a 7r-ring can 



2L\ $2 



LI C 

7T$n 



■ COS 



2tt$ 

$0 



(8) 



where $ is the magnetic flux included in the ring, and 
<fr e is the flux externally applied to the ring. At tempera- 
tures important in the cooling process, close to the Nb su- 
perconducting transition temperature T c , field screening 
effects can be neglected, since the magnetic penetration 
lengths are larger than the size of the sample. Under 
these conditions, we estimate the ring self- inductances 
(for 2.7 ^m junction width rings spaced by 11.5/im) to be 
«3.9pH, and the mutual inductance M between nearest 
neighbor rings to be «0.025pH. For a single junction n- 
ring spontaneous supercurrents flow when f3 — 2itLI c /&q 
is just greater than 1. The ring ordering process will oc- 
cur in the regime (3 — 1 << 1. In this limit, and in the 
limit $ a /$o << 1) where $ Q is the magnetic flux applied 
per ring, the potential barrier to flipping the sign of the 
ring supercurrents in the i th ring can be approximated 

by 



E u = 3(1 - l/#) 2 / c $ /47r, 



(9) 



with /3* = (3 + 47r<J> Q cri/$0j where m = ±1 is the spin of 
the i th ring, and 



where $ e is the externally applied flux and 
mutual inductance between the i th and j th rings, 
same limits, the value of the spontaneous flux <J> m 
ated by the i th ring is given by 



$ 



2vr 2 



(1 - l/#). 



(10) 

is the 
In the 
gener- 

(11) 



We take the junction critical current I c (0)=0.55mA for 
these rings, and the self-inductance of the rings is as- 
sumed to be temperature independent. A measure of the 
spin-spin coupling energy E c is the difference between 
the energy barrier for all nearest neighbor spins up, mi- 
nus that energy barrier if one of the neighbor spins is 
down. 

Fig. 1111 plots the calculated barrier to thermally acti- 
vated flipping Eb, nearest neighbor coupling energy E Cl 
and the thermal energy ksT as a function of temper- 
ature for the honeycomb ring array of Fig. Both 
Eft and E c increase as the temperature is lowered un- 
til the flipping freezes out at a temperature Tf such 
that E b /k B T f « - In^/ikB^Tf^dT/dt), where 7 = 
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of our simulations do not depend sensitively on whether 
periodic or non-periodic boundary conditions are used. 
For each iteration cycle, the probability of a spin flip Pj, 
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FIG. 12: Bond order as a function of reduced temperature 
T/T c calculated using Monte-Carlo techniques as described 
in the text, for various cooling rates rj. 



dEb(Tf)/dTm^ For the rings of Fig. El we estimate this 
occurs at (T max - T f )/T c = 3xl0~ 2 , with T max /T c 
1 — 1/Po and Pq — 2irLI c (0)/$o. At this temperature 
Ec/ksTf « 1.2. Therefore antifcrromagnetic coupling is 
energetically favored in isolated 7r-rings, but not nearly 
as strongly as in facetted junctions. 

As the temperature is lowered through T c , the junc- 
tion critical currents increase until two distinct circulat- 
ing states become allowed for (3 > 1 . Thermally activated 
switching between these states has a transition rate 



2ttt ( 



-E bi /k B T 



(12) 



where (3 C — 2ttI c R 2 C /<I>o is the junction hysteresis pa- 
rameter, and r = ^o/2ttI c R is the junction charac- 
teristic time^ We take the junction capacitance C = 
5.4 x lCr 13 F, J c (0) = 0.55mA, J C (T) = 7 C (Q)(1 - T/T c ), 
and R—lil for the 2D arrays presented in this paper. We 
have taken the limits where the d.c. shielding currents 
are much smaller than the critical currents, the spin- flip 
transition times are much shorter than the time between 
spin- flips, and f3 c « 1. We can isolate the tempera- 
ture dependence of the prefactor of Eq. ^] by writing 
(3 c /{2ttt) = 47r 2 / c 2 (0) 2 i? 3 C(l - T/T c ) 2 /$ 2 . It is of inter- 
est to note that the power dissipated per ring near the 
"freezout" temperature is of order ksT c /T ~ lO -11 ^. 

To model the cooling process in our 2D ring arrays 
we use a Metropolis Monte Carlo^ 3 . simulation: A 30x30 
element array is set up with the same geometry as the 
experimental arrays. Boundary effects are minimized by 
using periodic boundary conditions: rings at one edge 
of the array are treated as if their nearest neighbors are 
the corresponding rings at the opposite edge. The results 



Pi = (1 - T/T c f 



-E bi /k b T 



(13) 



is calculated for each ring. A random number between 
and 1 is generated. If this number is less than Pj, the spin 
of the ring is flipped. The process is repeated throughout 
the array, and the temperature is gradually reduced until 
no more spin flips occur. Figure H2I shows results from a 
Monte-Carlo simulation of the honeycomb lattice of Fig. 
13 using the ring parameters of Table I, plotting the bond 
order a as a function of the reduced temperature T/T c , 
for various values of 77, the change in T/T c per iteration 
cycle. Note that in this case the maximum temperature 
T m ax for which spontaneous magnetization is observed is 
T max /T c = 1 — l//?o = 0.846. As the temperature is re- 
duced the spin flipping probabilities decrease, and anti- 
ferromagnetic ordering gradually occurs. Furthermore 
note that the freezing temperature is well below the tem- 
perature at which E ~ k T. The effective cooling rate 
of the simulation is given by dT/dt = T]T c f3 c (0) /2ttt(0) = 
2.2 x 10 12 77 K-sec -1 . Therefore the slowest cooling rate 
simulated is about 1 x 10 4 K/sec, much faster than the 
experimental cooling rates. Although it would in princi- 
ple be possible to use sufficient computer time to match 
the modelled cooling rates to experiment, the modelling 
indicates stronger anti-ferromagnetic ordering than is ob- 
served, even at these very fast cooling rates. 

Fieure lTSl shows experimental results for the final bond 
order for the honeycomb array of Fig. [fJl for various cool- 
ing rates. Experimentally the effect of cooling rate on 
final bond-order is weak. Also shown in this Figure are 
the results of our Monte-Carlo simulations for the honey- 
comb array. The points labelled M nn = 25 fH use all of 
the parameters in the simulation as calculated following 
Table I. The points labelled M nn = 6fH used a reduced 
value of the ring-ring mutual inductance. Also shown 
are curves comparing the predictions if only the near- 
est neighbor spin-spin interactions are included, with re- 
sults if 2nd and 3rd nearest neighbors are also included. 
The mutual inductance between a ring and its further 
neighbors (e.g. 2nd nearest and 3rd nearest) is taken 
to be scaled by the cube of the ratio of the relevant 
center-to-center distances. Qualitatively, the inclusion 
of next nearest neighbors reduces the tendency towards 
anti-ferromagnetic ordering. However, it appears that 
the inclusion of 3rd nearest neighbors does not change 
the results much further. Similarly, reducing the strength 
of the spin-spin coupling by a factor of 4 also reduces the 
tendency to order. However, in all cases it appears that 
if the simulation were to be extended to sufficiently slow 
cooling rates to match experiment, long-range ordering 
would result. 

This failure of our arrays to order implies that there is 
a source of disorder that we haven't taken into account. 
There are three additional sources of disorder that we 
have considered. The first is a distribution in Nb critical 
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FIG. 13: Low temperature limit ol the bond order, as deter- 
mined experimentally with SQUID microscopy imaging (solid 
symbols) for the honeycomb lattice of Fig. |S| and as calcu- 
lated for the same lattices using Monte-Carlo techniques as 
described in the text (open symbols), for two values of the 
mutual inductance between nearest neighbor rings. The cal- 
culations labelled 1st include only nearest neighbor spin-spin 
coupling, those labelled lst+2nd also include second nearest 
neighbors, and those labelled lst+2nd+3rd also include third 
nearest neighbors. 
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temperatures, and the second is a distribution in junc- 
tion critical currents. In each case, our simulations in- 
dicate that an unrealistically large distribution width is 
required to significantly effect the predicted final a value 
vs. cooling rate curves. In all the cases we have ex- 
plored numerically: reducing Ef,, reducing E c , increasing 
the width of the T c distribution, or increasing the width 
of the I c distribution, the a vs. cooling rate curves can 
be shifted vertically in Fig. ^| but the slope of these 
plots does not change qualitatively. A final factor that 
might inhibit ordering of our arrays is if the cooling rate 
is not uniform. However, inspection of Fig. 1121 indicates 
that fast jumps in temperature of order 0.01T C i=s lOOmi-T 
would be required to significantly affect the final bond- 
orders. It seems unlikely that our temperature sweeps 
are that non-uniform. It is therefore difficult to see how 
the experimental data can be fit using the present simu- 
lations, and some factor that we haven't considered cor- 
rectly is causing the failure of our arrays to achieve long 
range order. 

Nevertheless, the experimentally determined values for 
the bond-order a can be modelled fairly well, if one uses 
the cooling rate as a fitting parameter. The left-hand 
panel of Fig. 114b summarizes cooling experiments for 
all of the 2D arrays of Fig. |SJ These data were taken 
after cooling the arrays at rates between lmK/sec and 
lOmK/sec. Although long-range anti-ferromagnetic or- 
der was not observed in these arrays, they did show 
strong anti- ferromagnetic correlations, as evidenced by 
the large negative bond orders. Note that there are no 
qualitative differences between the geometrically frus- 
trated and non-frustrated lattices in their tendency to 
anti-ferromagnetically order. For comparison, the 0-rings 
of Ref . |3i never showed bond-orders more negative than 
(7 = —0.2, whereas our 7r-ring arrays attained bond or- 
ders as negative as a = —0.3. The right hand panel 
shows the results of Monte-Carlo simulations for these 
arrays, using the ring parameters in Table I, except for 
Mnn = 6fH instead of M nn = 25fH, and with a cooling 
rate 1 x 10 6 K /sec. In both panels the ideal curve for an 
unfrustrated lattice (solid line) and a frustrated lattice 
(dashed line) are also indicated. The experimental re- 
sults can be qualitatively modelled assuming a spin-spin 
coupling constant about 4 times smaller than calculated, 
and a cooling rate 10 8 times faster than the experiment. 



FIG. 14: (a) Measured dependence of the bond order a upon 
fraction of spin- up rings x+, for the four different types of 
arrays shown in Fig. HJ obtained by cooling in various fields. 
The solid line is the theoretical maximum negative bond order 
for a geometrically unfrustrated array; the dashed curve is 
that for a frustrated array. Some of the experimental points 
at very low and high x+ fall below the theoretical maximum 
negative bond order because of edge effects due to the finite 
array sizes imaged, (b) Results of a Monte-Carlo simulation 
of the cooling process for these arrays, as described in the 
text. 



V. CONCLUSIONS 

We have shown that it is possible to use large arrays of 
photolithographically patterned 7r-rings as a model spin 
system. Half-fluxon Josephson vortices in electrically 
connected ID arrays, in the form of zigzag — 7r facetted 
junctions, order strongly anti-ferromagnetically through 
the superconducting order parameter phase. Electrically 
isolated ID and 2D arrays order much less strongly. The 
2D 7r-ring arrays show stronger anti-ferromagnetic cor- 
relations than reported previously for 0-ring arrays, but 
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not long-range order beyond a few lattice constants. We 
can understand the anti-ferromagnetic phase coupling of 
the ID zigzag junctions by solving the Sine-Gordon equa- 
tions. Although some features of the ordering in the 2D 
arrays can be understood using Monte-Carlo simulations, 
unrealistic parameters must be used, and we do not un- 
derstand why these arrays do not show long range order. 



VI. ACKNOWLEDGEMENTS 



We would like to thank D.H.A. Blank, R.H. Koch, K.A. 
Moler, G. Rijnders, and H. Rogalla for useful discussions. 
We also acknowledge the Dutch Foundations FOM and 
NWO for supporting this research. 



1 B.S. Deaver and W.M. Fairbank, Phys. Rev. Lett. 7, 43 
(1961). 

2 R. Doll and M. Nabauer, Phys. Rev. Lett. 7, 51 (1961). 

3 D. Davidovic et al., Phys. Rev. Lett. 76, 815 (1996). 

4 D. Davidovic et al., Phys. Rev. B 55, 6518 (1997). 

5 L.N. Bulaevskii, V.V. Kuzii, and A. A. Sobyanin, JETP 
Lett. 25, 290 (1977). 

6 V.B. Geshkenbein and A.I. Larkin, JETP Lett. 43, 395 
(1986). 

7 V.B. Geshkenbein, A.I. Larkin, and A. Barone, Phys. Rev. 
B 36, 235 (1987). 

8 DA. Wollman, D.J. Van Harlingen, W.C. Lee, D.M. Gins- 
berg, and A.J. Leggett, Phys. Rev. Lett. 71, 2134 (1993). 

9 DA. Brawner and H.R. Ott, Phys. Rev. B 50, R6530 
(1994). 

10 C.C Tsuei et al., Phys. Rev. Lett. 73, 593 (1994). 

11 A. Mathai, Y. Gim, R.G Black, A. Amar, and F.C. Well- 
stood, Phys. Rev. Lett. 74, 4523 (1994). 

12 H.J.H. Smilde et al., Phys. Rev. Lett. 88, 057004 (2002). 

13 H.J.H. Smilde, H. Hilgenkamp, G. Rijnders, H. Rogalla, 
and D.H.A. Blank, Appl. Phys. Lett. 80, 4579 (2002). 

14 A.V. Andreev, A.I. Buzdin, and R.M. Osgood, III, Phys. 
Rev. B 43, 10124 (1991). 

15 V.V. Ryazanov, VA. Oboznov, A.V. Veretennikov, and A. 
Yu. Rusanov, Phys. Rev. B 65, 020501(R) (2001). 

16 A. Bauer et al., Phys. Rev. Lett. 92, 217001 (2004). 

17 H. Hilgenkamp et al., Nature 422, 50 (2003). 

18 F.P. Rogers, Master's Thesis, MIT, Boston, (1983). 



19 L.N. Vu and D.J. van Harlingen, IEEE Trans. Appl. Su- 
percond. 3, 1918 (1993). 

20 R.G Black et al., Appl. Phys. Lett. 62, 2128 (1993). 

21 J.R. Kirtley et al., Appl. Phys. Lett. 66, 1138 (1995). 

22 D. Xu, S.K. Yip, and JA. Sauls, Phys. Rev. B 51, 16233 
(1995). 

23 H. Susanto, SA. van Gils, T.P.P. Visser, Ariando, H.J.H. 
Smilde, H. Hilgenkamp, Phys. Rev. B 68, 104501 (2003). 

24 V.G. Kogan, J.R. Clem, and J.R. Kirtley, Phys. Rev. B 
61, 9122 (2000). 

25 J.R. Kirtley, K.A. Moler, and D.J. Scalapino, Phys. Rev. 
B 56, 886 (1997). 

26 E. Goldobin, D. Koelle, and R. Kleiner, Phys. Rev. B 66, 
100508(R) (2002). 

27 E. Goldobin, D. Koelle, and R. Kleiner, Phys. Rev. B 67, 
224515 (2003). 

28 A. Zenchuk and E. Goldobin, Phys. Rev. B 69, 024515 
(2004). 

29 E. Goldobin, N. Stefanakis, D. Koelle, and R. Kleiner, 
Phys. Rev. B 70, 094520 (2004). 

30 H.J.H. Smilde, Ariando, D.H.A. Blank, H. Hilgenkamp, 
and H. Rogalla, Phys. Rev. B 70, 024519 (2004). 

31 J.R. Kirtley, GG Tsuei, and F. Tafuri, Phys. Rev. Lett. 
90, 257001 (2003). 

32 CD. Tesche, J. Low. Temp. Phys. 44, 119 (1981). 

33 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. 
H. Teller, and E. Teller, J. Chem. Phys. 21, 108 (1953). 



